Quantum Critical Dynamics Simulation of Dirty Boson Systems 
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Recently the scaling result z = d for the dynamic critical exponent at the Bose glass to superfluid 
quantum phase transition has been questioned both on theoretical and numerical grounds. This 
motivates a careful evaluation of the critical exponents in order to determine the actual value of z. 
We study a model of quantum bosons at T — with disorder in 2D using highly effective worm 
Monte Carlo simulations. Our data analysis is based on a finite size scaling approach to determine 
the scaling of the quantum correlation time from simulation data for boson world lines. The resulting 
critical exponents are z = 1.8 ± 0.05, v — 1.15 ± 0.03, and r\ = —0.3 ± 0.1, hence suggesting that 
z = 2 is not satisfied. 
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Quantum phase transitions (QPT) occur at zero tem- 
perature and produce new and important physics com- 
pared to "classical" phase transitions at finite tempera- 
ture [H, 0. In particular, presence of quenched disorder 
leads to new universality classes without direct classi- 
cal counterparts, where lack of space-time symmetry can 
lead to nontrivial scaling properties. Such phenomena are 
of great current interest and present considerable theo- 
retical and experimental challenges P, Q ■ 

A prototype QPT with disorder is the 2D boson su- 
perfluid to insulator transition in the presence of ran- 
dom substrate disorder. The disorder localized insulating 
phase is a gapless phase called the Bose glass. This tran- 
sition is relevant for experiments on ultrathin granular 
superconducting films, Josephson junction arrays, super- 
fluid helium films, and cold bosons in optical lattices with 
disorder [1, 0] . A remarkable result of the theory is the 
relation z — d, where d is the number of spatial dimen- 
sions 3]. This scaling result was believed to be exact, 
but has been questioned recently both analytically and 
numerically [5H7|]. The result z — d is derived by requir- 
ing the contribution to the compressibility k from the 
singular part of the free energy to be a nonsingular func- 
tion across the transition. However, if k instead comes 
from the analytic part of the free energy no restriction 
on z follows and the relation z = d does not have to 
hold @. In ID z = 1 is fulfilled 0, but Ref. [6] finds 
that this is unrelated to the mechanism that keeps the 
compressibility finite through the Bose glass-superfluid 
transition. 

The task of determining the quantum dynamical expo- 
nent at the disordered boson QPT to test the validity of 
the relation z = 2 in 2D has been studied previously. An 
often used approach has been to assume the value z = 2 
and then test if scaling can be obtained by fitting other 
parameters to numerical data. This approach produces 
seemingly good scaling results for the system sizes tested 
[§-11 1, but does not rule out that a calculation without 
a priori assumptions might give a different result. A re- 
cent simulation study reports z » 1.4 0, but this result 



mig ht be affected by the limited disorder averaging used 
ll| . Renormalization approaches have also been used to 
determine z [12], but have not yet settled Q. Thus the 
validity of the result z — d is unclear and further tests 
are required. 

In this paper we perform large scale Monte Carlo (MC) 
simulations to determine z and other critical exponents 
at the Bose glass transition of the dirty boson model in 
d — 2 dimensions. We extend previous simulation results 
in several ways. We use extensive disorder averaging, and 
larger system sizes than in most previous studies, which 
turns out to be crucial. A highly effective worm algo- 
rithm is used that permits efficient averaging over con- 
figurations with different boson winding numbers 13[ . In 
order to locate the QCP and study dynamical scaling, a 
suitable function of the winding number is constructed 
that has a maximum value when the system size in the 
time direction is proportional to the correlation time. Fi- 
nite size scaling of the maximum gives a direct route to 
calculating z and other critical exponents, without any a 
priori assumptions on z. The results display significant 
corrections to scaling for small system sizes that com- 
plicates determination of the exponents. Our estimate, 
z = 1.8 ± 0.05, suggests that the dynamic exponent is 
smaller than given by the relation z = d for d = 2. 

Dirty boson model - The imaginary-time path-integral 
representation of the 2D Boson Hubbard model with 
nearest neighbor hopping, on-site charging energy, and 
a disordered chemical potential can be mapped to a link- 
current model convenient for simulation [9|. The link- 
current model assumes only phase fluctuations of the or- 
der parameter and neglects amplitude fluctuations, has 
isotropic space-time couplings, and uses the Villain form 
of the potential Q . Such details are not expected to al- 
ter the universality class of the QPT. The Hamiltonian 
of the link-current model is 
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Here i = (r, t) denotes the sites of a (2 + l)-dimensional 
simple cubic spacetime lattice of size L x L x L T with 
periodic boundary conditions in space and time direc- 
tions, and 5 = x,y,T denotes the coordinate directions. 
The integer link variables jf represent boson current 
variables on the links extending from the site i in the 
(^-direction. The variables are subject to the divergence- 
free constraint V- J = 0, which means that the worldlincs 
have no open ends. K is a coupling constant. Disorder 
is modeled as a quenched on-site potential which is ran- 
dom in space but constant in the time direction, with a 
uniform distribution in \v r \ < 1/2. The chemical poten- 
tial is here fixed to [i = 1/2 which means half filling of 
bosons on average. The transition of this model repre- 
sents the generic universality class of the disorder driven 
boson localization QPT. 

Next we introduce the two main quantities of interest 
in our simulations. The mean square winding number is 
defined as 
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The bracket (• • •) indicates average with respect to J- 
current configurations, and [• • • ] indicates the quenched 
disorder average. The spatial mean square winding num- 
ber measures fluctuations in the number of times the 
worldlines wind across the sample, and is proportional 
to the superfluid density [1J]. It can thus be used to de- 
tect the boson superfluid to insulator QCP. The tempo- 
ral winding number fluctuations (including subtraction 
of the average boson number) correspond to the boson 
compressibility k [9j. The gapless nature of the Bose 
glass produces a smooth nonzero compressibility across 
the transition [3J. From now on we will only consider 
spatial winding number fluctuations, and form W 2 = 
(W 2 + W 2 )/2. The Greens function G(r -r',r- r') = 



')} can be used to define the uniform order 



parameter susceptibility \ = G(k = 0, w = 0) @, 

Monte Carlo simulations - Our MC simulations use 



the classical lattice worm algorithm 13]. For each dis- 
order realization the simulation was started in the J- 
current configuration that minimizes H in Eq. (p}. The 
simulations used more than 1500 MC sweeps to reach 
equilibrium, followed by equally many sweeps to collect 
data for the averages. Here a MC sweep is defined as 
3L 2 L T link variable update attempts. Measurements are 
taken every time the worm closes. The winding number is 
given by the number of times the world lines wrap around 
the sample, and the susceptibility is the average number 
of update attempts per closed loop configuration [Icj . 
We tested for equilibration by monitoring disorder aver- 
ages of the winding number fluctuations and of the sus- 
ceptibility calculated using different numbers of warmup 
sweeps. An example is shown in the inset in Fig. [TJ The 
results become independent of the initial configuration 



after about 500 warmup sweeps. The quenched disorder 
averaging used between 10 4 — 10 5 samples of the ran- 
dom potential, where more disorder averaging was used 
around the critical point. Statistical error bars on the 
data points were estimated by fluctuations in the disor- 
der averages. 

Finite-size scaling methods - The basic scaling as- 
sumption is that the correlation length and time di- 
verge at the transition as £ ~ \k\~ v and r ~ £ z , where 
k = (K — K c )/K c , K c is the critical coupling, v is the 
correlation length exponent, and z is the dynamic expo- 
nent. The winding number fluctuation is dimensionless 
and therefore scale invariant at the transition. We as- 
sume the following finite size scaling (FSS) ansatz for 
the winding number fluctuation 



W 2 (K, L, L T ) = W 2 {L l/u k, a T 



and for the susceptibility 

x(K,L,L T ) = L 2 -r>x(L 1 ^k,a T ) 
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where W 2 and x are scaling functions, and a T = L T /L Z 
is the aspect ratio. The aim is to estimate the critical 
exponents z,v,rj and scaling functions by fitting these 
expressions to numerical MC data for finite L, L T . 

FSS analysis greatly simplifies if the scaling functions 
can be reduced to functions of only one variable by taking 
the other variable to be constant. Taking the first vari- 
able L x l v k to be constant means keeping K = K c , which 
is a priori unknown, while keeping the second variable 
constant requires knowledge of z. Most previous stud- 
ies have therefore assumed the value z — 2 and selected 
system sizes for simulations given by L T = const x L? . 
Clearly this approach is not available if the value of z is 
unknown. 

The idea is now to, without assuming knowledge of K c 
and z, construct a characteristic scale L* for each given 
K, L, which scales as L* ~ r ~ L z for K = K Cl where 
r is the correlation time. The winding number fluctua- 
tion is a monotonically increasing function of L T for fixed 
K, L. For L T 3> r the worldline fluctuations separated 
by times greater than the correlation time r decorrelate, 
and then the winding number fluctuation must increase 
linearly with L T . Thus the quantity W 2 / L T approaches 
a constant value for L T 3> t. Dividing once more gives 
W 2 /L 2 , which has a maximum at a characteristic L*, 
and goes to zero for L T 3> t, where the star indicates 
the value at the maximum. We will find these maxima 
very useful in the scaling analysis [HI. A convenient 
scaling form is produced by replacing L T in W 2 /L 2 by 
a T = L T /L Z . We thus introduce 



$(K,L,L T 
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This FSS relation is used below to estimate the critical 
coupling K c and the exponents z, v. We verified that our 
approach reproduces known exponents for pure models. 
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FIG. 1. Selection of Monte Carlo results for the winding num- 
ber fluctuation divided by L 2 as a function of L T . Solid curves 
are polynomial fits to the data curves, from which the loca- 
tions L* and sizes (W 2 /L 2 )* of the maxima can be deter- 
mined. Inset: Equilibration test for L — 40, L T = 240, K — 
0.2477. 
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FIG. 2. Maximum values $* = (W 2 /a 2 -)* vs. system size 
L for different couplings K. For z — 1.78 the data obeys 
$* = const at K = K c for L > 16, indicated by the horizontal 
dashed line, which estimates K c — 0.2477. The lower dashed 
line corresponds to z — 2, which approximately describes the 
data for small sizes L < 16 at K = 0.246. 



Results - First we locate the critical coupling K c and 
the dynamic exponent z by FSS analysis of MC data for 
the winding number. Figure [1] shows examples of maxima 
of the quantity W 2 /L 2 T . The amplitude (W 2 /L 2 )* and 
location L* of the maxima can be straightforwardly com- 
puted by polynomial fits to the MC data curves. Better 
accuracy is obtained in the estimates for (W 2 /L 2 )* than 
for L*. The maximum values scale as (W 2 /L 2 )* ~ L~ 2z 
at K = K c . However it is more convenient to plot the 
quantity $* = {W 2 /a 2 T )* of Eq. ©, and look for the 
scaling <f>* ~ L° at K — K c , which is shown in a log-log 
plot in Fig. [2j In the fi gure the of value z enters through 
a* = L*/L z , and has been adjusted to make $* = const 
at K = K c for system sizes L > 16, marked with the hor- 
izontal dashed line. This produces the estimates z re 1.8 
and K c re 0.2477. For K ^ K c the data curves clearly 
splay out, away from a critical power law. For L < 16 
deviation from power law behavior is obtained, which in- 
dicates the presence of corrections to scaling in these data 
points. In Fig. [2] we also note that the choice z = 2 gives 
an approximate description of the data for K = 0.246 
for small system sizes, L < 16, which is indicated by the 
lower dashed line, in agreement with Ref. [l(| . As a con- 
sistency test, a similar analysis was done for the location 
L* of the maxima using the relation L* ~ L z at K = K Cl 
leading to similar results. 

Figure [3] A shows the maxima of the function $* of Eq. 
([5]), with a = L T / L z for z re 1.8. The data curves for L > 
16 intersect at K c , but for smaller sizes scaling deviations 
are present, and these will be further discussed below. 
The correlation length exponent is readily estimated by 
computing the derivatives d$*/dK\K=K c ~ L 1 '", and a 
polynomial fit to the MC data gives v re 1.15. The FSS 
data collapse produced by using this value for v is shown 




FIG. 3. A: Intersection plot for the scaled winding number 
function maxima $*, showing an intersection at K c = 0.2477. 
B: FSS data collapse of the data in A obtained for v = 
1.15,L > 16. C: Intersection plot for the scaled susceptibility 
x/i 2- " 7 evaluated at a T = 0.35. The data curves for large 
sizes roughly intersect at K — 0.2477 for r\ — —0.29, but with 
much larger corrections visible for small system sizes than for 
$*. D: FSS data collapse of x/i 2 "" for v = 1.15, L > 16. 

in Fig. OB for L > 16. 

To estimate the correlation function exponent r\ we use 
the susceptibility x given by Eq. . We fix the aspect 
ratio to a T — 0.35 which correspond to the maxima 
at criticality. The value of \ at this aspect ratio was de- 
termined by a polynomial fit to nearby MC data. From 
X ~ L 2 ^ v we estimate r\ re —0.3 for L > 16. Figure [3] C 
shows a corresponding intersection plot for the quantity 
x/L 2 ~ v , which becomes size independent at K = K c ac- 
cording to Eq. Q. A FSS collapse assuming v = 1.15 is 
shown in Fig.[3]D. Note that the deviations from scaling 
for small system sizes in Fig. [3] C are substantial, and 
hence the uncertainty in the estimate of r\ is consider- 
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FIG. 4. FSS data collapses of MC data for the scaled winding 
number fluctuation $ = W 2 /a 2 and susceptibility x/^ 2 ^ 
as functions of a T = L T /L Z for K c = 0.2477, z = 1.8, 77 = 
—0.29. Solid curves are polynomial fits to the data. Inset: 
Dependence of z on the range of system sizes used in the 
estimate. L m i n indicates the smallest in a sequence of size 
quadruplets in L = 8, 10, 12, 16, 20, 30, 40, 60 used to estimate 
z, except for L min = 30 which indicates sizes 30, 40, 60. 



able. The scatter among the intersection points can be 
reduced by assuming a scaling correction proportional to 
L~ u with u « 1, but the accuracy of the data is insuffi- 
cient for detailed estimates. 

Finally we systematically study the system size depen- 
dence of the estimated exponents and estimate errors. 
This final calculation does not involve the maxima, and 
thus avoids any errors in their determination. A double 
polynomial expansion is done of the scaling functions in 
Eq. (JSJ) in both arguments. The parameters are deter- 
mined by ^-minimization of the RMS deviations of the 
MC data points from the scaling function. We performed 
several fits for MC data points selected from different in- 
tervals in the range 0.2 < a T < 1.2 in order to verify 
the stability of the results. To study system size trends 
of the results, fits were made for a sequence of system 
size quadruplets in L = 8,10,12,16,20,30,40,60. The 
result for z is shown in the inset of Fig. 01 The dis- 
played trend agrees with the one indicated in Fig. [2j For 
the fit with L = 16, 20, 30, 40, 0.2 < a T < 0.5, we get 
X 2 /DOF « 0.8. Our final estimates including error esti- 
mates based on statistical errors determined by the boot- 
strap method combined with average variations from the 
dependence on the a T -interval included in the fits are 
K c = 0.2477 ± 0.0002, z = 1.8 ± 0.05, v = 1.15 ± 0.03, 
and r] = —0.3 ± 0.1. Other critical exponents can be 
estimated from these values using scaling laws. 

Discussion - Analysis of our MC data of the 2D bo- 
son localization transition by disorder revises previous 
estimates of the critical exponents. In particular the dy- 



namic critical exponent is estimated to z = 1.8 ± 0.05, 
which suggests that z = d is not fulfilled in d = 2, al- 
though the values are close. Our results clarify how most 
previous simulations appear consistent with z — 2. For 
small system sizes z — 2 works quite well, but including 
larger sizes reveals corrections to scaling making z = 1.8 
a better estimate. Our estimates are quite different from 
those of Ref. 0] , which we believe may be explained by 
their smaller disorder averaging and uncertainty in their 
location of the quantum critical point. The prediction 
of a universal conductivity at the transition is indepen- 
dent of the value of z @. However, actual estimates of 
the universal value of the conductivity indirectly depend 
on the value of z, and should be reexamined in the light 
of the present results. A better analytic understanding 
of the quantum critical dynamics as well as further ex- 
perimental measurements probing these issues would be 
welcome. 
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